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A simple scenario of the formation of geological landscapes is suggested and the respective lattice 
model is derived. Numerical analysis shows that the arising non-Gaussian surfaces are characterized 
by the scale-dependent Hurst exponent, which varies from 0.7 to 1, in agreement with experimental 
data. 
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Rough interfaces around us, such as Earth's surface 
jl], 0], surfaces of deposited films j|, wetting fronts jjj, 
cloud perimeters fracture surfaces ||, etc., are com- 
mon objects, the properties and formation of which have 
been studied for several decades. In some cases, signif- 
icant advances in theoretical understanding have been 
achieved. In particular, this applies to the surface growth 
processes, the analysis of which has led to a wide variety 
of kinetic roughening models, cf. [Q, |). However, many 
processes leading to rough surfaces, eg. the formation of 
fractures and Earth's landscapes are less understood. 

The formation of the Earth's surface is a complex pro- 
cess, affected by various phenomena, such as seismic and 
tectonic activity, erosion, sedimentation, etc. These phe- 
nomena, in their turn, can be of diverse nature. Thus, 
erosion can be caused by meandering rivers, oceanic and 
atmospheric influence, by the motion of ice, avalanches, 
and so on. Furthermore, the physical properties of the 
ground vary in a very wide range. Incorporating all this 
diversity into a concise embraceable mathematical model 
is a hopeless task. However, the scale-invariant proper- 
ties of the geologic landscapes appear to be surprisingly 
universal: in a reasonable approximation, they are typi- 
cally self-affine, with the Hurst exponent ranging between 
H s» 0.7 and 0.9 0, g]. Therefore, it is natural to expect 
that there is a simple, universal and robust mechanism 
leading to such surfaces. In what follows we show that 
such a mechanism can be provided by the competition 
of erosion and tectonic activity: the model of gradient- 
limited surfaces incorporates these two effects in their 
simplest form and leads to realistic landscapes. 

Most models of geological landscapes are based on the 
evolution of river networks §, y, [l|. The evolu- 
tion of rivers plays undoubtedly an important role in the 
formation of landscapes, but is not able to increase the 
height of the mountains and will not lead to self-affinc 
surfaces. The only attempt of constructing a robust self- 
affine model of Earth's surface has been made by Mandel- 
brot g, §. He modelled roughening due to tectonic ac- 
tivity, and his method can be outlined as follows. Inside 
a polygon (Earth's surface), a random point is coined. 
Through that point, a line of random direction is drawn. 
This "fault line" divides the polygon into two parts, one 
of which is elevated (with respect to the other) by a unit 
height. The procedure is repeated N — ► oo times. The 
Brownian growth of height differences is eliminated by 



normalizing the surface height to ^f~N. This results in 
a self-affine surface with H — 0.5. In order to address 
the discrepancy between the model and empirical values 
H « 0.7 - 0.9, the model has been generalized by re- 
placing the Heaviside profile of the "fault" by a profile 
with singularity (so that the height change is given by 
Ah = assigns, where —0.5 < a < 0.5 and the x-axis is 
perpendicular to the "fault line"). 

The tectonic activity and formation of faults, as cap- 
tured by the Mandelbrot's model, plays certainly an 
important role in the evolution of the Earth's surface. 
Meanwhile, the singular shape of the fault profile is ar- 
tificial, with no physical motivation. Besides, there are 
no physical processes which would normalize the surface 
height to the number of faults. Instead of that, the ba- 
sic effect reducing the height differences is erosion. As 
mentioned above, erosion itself is a very complex phe- 
nomenon which, in particular, has been addressed by the 
models of the evolution of river networks [g [j^, [TJ, fTJ . 
The excessively detailed erosion models, however, are not 
suited for revealing the most generic aspects of landscape 
roughening. Therefore, we opt for the simplest possible 
approach and assume that effectively, erosion imposes an 
upper limit to the modulus of the gradient of the sur- 
face height. More specifically, we assume that, as soon 
as a slope becomes steeper than a threshold value, the 
excess of the height drop is spread over the neighboring 
regions. Such a smoothing of too steep slopes can be 
accomplished, for instance, by avalanches. 

To begin with, let us define the model of gradient- 
limited surfaces for a continuous medium. A random 
point P and direction t define a "fault line" inside a 
polygon of diameter L 1 . This line divides the polygon 
into two parts, one of which (leftmost, with respect to the 
direction t) is elevated by a unit height. The height drop 
is spread over the nearest neighboring regions in such 
a way that the modulus of the local gradient remains 
everywhere below a threshold value, i.e. |VV>| < 1. If 
the fault line goes through a region of a saturated slope, 
the avalanches can affect large areas. Then, the actual 
elevation (the change of surface slope) will take place 
far from the fault, at the edges of the saturated slope. 
These edges will be referred to as the elevation lines. The 
procedure is repeated N — ► oo times. Note that unlike 
the Mandelbrot's model, this model has a lower cut-off 
scale (the ratio of the height drop at a single fault line 



and threshold gradient). 

The model can also be formulated on a lattice. A nat- 
ural basis for the lattice formulation is given by the six- 
vertex (restricted solid-on-solid) model || [f|. On the 
square lattice of the six vertex model, all the edges are 
marked with arrows so that each site has equal number of 
incoming and outgoing arrows, see Fig. 1. The arrows de- 
fine an incompressible flow, the streamfunction of which 
is our surface. Thus, each arrow represents a unit jump 
in the surface height. The mean slope is defined by the 
mean density of counter-directed arrows, which is lim- 
ited by the step of the grid; hence, the slope steepness is 
constrained automatically. 

Consider an oriented chain of arrows dividing the lat- 
tice into two parts (see Fig. 1). Swapping the direction 
of all the arrows of the chain is legitimate, because for 
all the affected sites, the number of incoming arrows is 
conserved. It corresponds to the lowering of one part of 
the surface with respect to the other part by two units. 
Therefore, directed chains of arrows can play the role of 
elevation lines. Gradient-limited surfaces are obtained 
as follows. A random site of the lattice P and a random 
direction r are coined, they define an aim line, the site of 
the "fault" . That part of the surface, which is leftwards 
to the aim line, is to be elevated. The elevation is ac- 
complished along such a directed chain of arrows, which 
follows as closely as possible the aim line. Similarly to 
the continuous case, if the aim line goes through a region 
of uni-directional arrows (saturated slope), the elevation 
line is forced to go around those regions. There are dif- 
ferent technical options, how to minimize the distance 
between the elevation line and aim line. In particular, 
the distance can be optimized locally and globally. How- 
ever, the scaling properties of the resulting surfaces are 
insensitive with respect to the particular choice. For our 
main series of simulations, we used a local algorithm. 
The elevation line was traced step-by-step, starting from 
the origin P. At each step, there are two possibilities 
to continue the line, because there are two outgoing ar- 
rows from each site (however, if the site has been al- 




FIG. 1: The algorithm of finding the elevation line (depicted 
by bold line). The aim line (dashed) is a random line through 
a randomly chosen (white) vertex. The elevation line is such a 
directed chain of arrows, which follows the aim line as closely 
as possible. Darker areas (squares) are lower. 
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FIG. 2: A gradient-limited surface, polygon size I/ max = 2049. 
Darker areas correspond to lower regions of the surface, black 
and white lines depict equi-distant level lines. Black line sur- 
rounded by white is an elevation line. 



ready visited, one of the outgoing arrows is occupied, 
and there is only one possibility left). That option is 
to be selected, which leads closer to the aim line. The 
(signed) departure from the aim line is easily tracked as 
s = Ax sin a — Ay cos a, where Ax and Ay are displace- 
ments along x and y axes, and a — the angle between the 
x-axis and the aim line. The elevation line is terminated 
as soon as it reaches the boundary of the polygon. The 
gradient-limited surfaces are obtained at the long-time 
limit, when the number of elevations exceeds the relax- 
ation time (which scales as the number of sites in the 
polygon, see below), and the initial shape of the surface 
becomes irrelevant. 

A simulation result is presented in Fig. 2. Observe 
the shape of the elevation line: in the region of a sat- 
urated slope, it has to depart far from the straight aim 
line. These regions are responsible for the long-range cor- 
relation of the surface height increments. Qualitatively, 
the saturated slopes are caused by accumulated excess 
of the "faults" of a certain direction. If there were no 
avalanches, that excess would fluctuate as the square root 
of the number of "faults" , tending to infinity. Therefore, 
the presence of large saturated slopes should not be sur- 
prising. As it will be shown below, in one-dimensional 
(ID) case, the accumulation phenomenon gives rise to 
H = 1 (i.e. typically, a saturated slope occupies the 
whole polygon). In 2D geometry, the accumulation ef- 
fect is weaker and leads to a more interesting scaling 
behavior. 

The one-dimensional (ID) version of the model is most 
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FIG. 3: Gradient-limited surfaces: numerical results. The 
numbers 1025, 513, 259, 127, and 65 indicate the edge length 
imax of the polygon. The finite difference approximation of 
the differential roughness exponent h is plotted versus the 
logarithmic relative scale A. There is no strict self-affinity 
of the surface. However, at the limit L max — > oo, there is 
an asymptotic dependance /i(A,I/ max ) — > h(X) (dotted line, 
dashed part is extrapolation). 

conveniently formulated as a spin exchange problem, and 
admits analytical solution. This analytic approach helps 
us to understand the features of the 2D-model. Sup- 
pose there is a sequence of spins, tpi = ±1. It is conve- 
nient to consider infinite periodic sequence, ifii+N = fit 
where i G Z and N is the period. The spins tpi can 
be interpreted as the increments of a self-affine curve 
ipj = 5Zi=o Vi ■ A random point k and a random spin 
increment v = ±2 are coined. The increment is to be 
added to the A:-th spin, or, if it is not possible (resulting in 
ipk = ±3), to the nearest suitable spin (i.e. to tpi = — is±l 
with minimal value of \l — k\). If there is no suitable spin 
at all, the next pair of k and v are coined. The proce- 
dure is repeated ad infinitum. Let us denote the relative 
number of positive spins by £. Then, the height drop 
of the above defined self-affine curve ipj at distance N 
is N 1 2£ — 1|. The quantity £ performs Brownian fluctu- 
ations, because at each time step, it is randomly incre- 
mented by ±A _1 . At the limit N — > oo, the probability 
density function n(£, t) evolves according to the diffusion 
equation, n t = Dn^ with D = (2N 2 )^ 1 and no flux at 
the boundaries, n^(0,t) = nj(l,t) = 0. The stationary 
solution n(£, t) = 1 allows us to calculate the delta vari- 
ance ((ip i+N - ^) 2 ) = N 2 r X (2^ - l) 2 ^ = N 2 /3, which 
corresponds to H = 1. The relaxation time of the spin 
exchange problem can be found as the diffusion time, 
r iV 2 (time is measured in the number of spin ex- 
changes) . The relaxation time of the 2D gradient- limited 
surfaces can be estimated in the same way, because the 
height difference between left and right edges of the poly- 
gon performs also nearly-Brownian fluctuations. 

For two-dimensional geometry, the simulations indi- 
cate that the gradient-limited surfaces are not strictly 
speaking self-affine. However, the data collapse of the 



simulation results can be achieved, when we introduce 
the scale-dependent differential Hurst exponent, defined 
as 

h(X) = -dlog(oi) /dlogL, A = logi/logL max . (1) 

Here cll is the height of the surface at the dis- 
tance L from the center of the polygon of size L max - 
The angular braces denote averaging over different re- 
alizations of the surface. In Fig. 3, the differen- 
tial Hurst exponent is approximated by h(X, £ max ) = 
log((a 2 ) / (a 2 +1 ))[log(L J /L l+ i)] _1 , where i and i + 1 are 
neighboring data-points and A = log(LiLi+i)/logL 2 nax . 
At the limit L max — » oo, the curves converge to 
the asymptotic function h(X). Note that at the ex- 
treme right-hand-side of the plot, the convergence of the 
/i(A, L max )-curves is not as good as elsewhere; this is 
explained by finite-size effects and by the fact that for 
A ~ 1, the finite differences fail providing an acceptable 
approximation for the derivative in Eq. (1). The rapid 
fall-off of the curves at the limit A — > 1 has the same ori- 
gin, and therefore, the values h < 0.65 are not reliable. 
For large scales (which are most interesting in the context 
of the Earth's surface) with A > 0.9, a better approach 
is to extrapolate the asymptotic curve, see dashed line in 
Fig. 2. The conclusion h « 0.7-0.9 for 0.5 < A < 1 is in 
a good agreement with the values recorded for geologi- 
cal landscapes, cf. Intriguingly, the roughness expo- 
nents of the fracture surfaces vary in the same range, cf 
i 0- 

Finally, the scaling law of the overall (edge-to-edge) 
height drop of the gradient-limited surfaces is given by 
the integral Hurst exponent H = J Q hdX: (a(L max ) 2 ) oc 
I/max! according to the simulations, H = 0.91 ± 0.01. 
This exponent equals to the area under the asymptotic 




FIG. 4: Differential fractal dimension of the contour loops 
of gradient-limited surfaces: numerical results. The numbers 
65-1025 indicate the polygon size Z/ max - The finite difference 
approximation d is plotted versus the scale A. Dotted line de- 
picts the asymptotic dependance d(A,L max ) — > d(X). Dashed 
line is calculated using the dependence h(X) (see Fig. 3), and 
the fractal dimension for Gaussian self-affine surfaces. 
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FIG. 5: The finite difference approximation k of the differ- 
ential scaling exponent characterizing the size-distribution of 
contour loops is plotted versus the scale A. Positive values of 
k indicate that there is an anomalously small number of small 
contour loops (as compared with the Gaussian surfaces). 

curve in Fig. 3, and this condition has been used to test 
the extrapolated curve in Fig. 3. 

We have observed that the differential Hurst exponent 
is scale-dependent. Such a generalized scale-invariance 
when critical exponents depend on scale, is not unique. 
For instance, similar behavior has been observed for cer- 
tain forest fire models In our case, the Hurst ex- 
ponent increasing towards small scales is caused by the 
presence of large areas of saturated slope. Indeed, con- 
sider a random pair of points. If the distance L between 
them is small, the points are likely to reside inside a sin- 
gle region of saturated slope. Hence, their average height 
difference scales almost as L implying fowl. On the 
other hand, larger saturated slopes are more rare than 
the smaller ones. Therefore, for a more distant pair of 
points, falling inside a single saturated slope is a rare 
event, and the conclusion h « 1 is no more valid. 

For Gaussian self-afhne surfaces, all the scaling ex- 
ponents of statistical topography are functions of the 
Hurst exponent H. However, the gradient-limited sur- 



faces are not Gaussian, as evidenced by the presence of 
large saturated slopes. Therefore, the differential Hurst 
exponent h(X) alone does not provide a complete de- 
scription of the surface. First we consider the differen- 
tial fractal dimension of the contour loops ("coastlines"), 
d = dlog(^) /dlogL. Here I is the length of a contour 
loop, and L its diameter. The numerical results are given 
in Fig. 4. For Gaussian surfaces, the fractal dimension of 
contour loops D(H) « 1.5 - 0.5H ff,|l7)- Dotted line in 
Fig. 4 is the asymptotic (L max — > oo) dependence d(X), 
and dashed line is the curve, calculated on the basis of 
the functions h(X) (from Fig. 3), and D(H). Evidently, 
D{h{X))^d{X). 

Even more pronounced mismatch between the Gaus- 
sian surfaces and the gradient-limited surfaces is observed 
for size-distribution of the contour loops. Let p(L) de- 
note the probability that a randomly chosen point be- 
longs to such a contour loop, the diameter of which is 
larger than L, but smaller than 2L. Then, for Gaus- 
sian surfaces we would expect that p(L) = L k , where 
k = D(H)- (2- H) < §, For gradient-limited sur- 
faces, the curves k(X, i max ) converge again to the asymp- 
totic dependence fc(A). This convergence with k(X) > 
is clearly observed up to the scale A « 0.85; above of that 
scale, convergence is very slow due to finite-size effects. 
The inequality k > means that, as compared with the 
Gaussian surfaces, there is significantly lesser number of 
small contour loops, which is quite a natural observa- 
tion. Indeed, the saturated slopes can be embraced by 
large contour loops, but leave almost no room for small 
ones. 

In conclusion, the new model of gradient-limited sur- 
faces leads to non-Gaussian surfaces which are charac- 
terized by the scale-dependent differential Hurst expo- 
nent. The latter varies from h « 1 for small scales, up to 
h w 0.7 for large scales. This is in a reasonable agreement 
with the experimentally observed roughness of real goo- 
logical landscapes. The relevance of the model to fracture 
surfaces deserves further analysis. 
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